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Optimal assignment of classes to classrooms [J , design of DNA microarrays [2], cross species 
gene analysis [3|, creation of hospital layouts [J|, and assignment of components to locations on 
circuit boards [a| are a few of the many problems which have been formulated as a quadratic 
assignment problem (QAP). Originally formulated in 1957, the QAP is one of the most difficult of 
all combinatorial optimization problems. Here, we use statistical mechanical methods to study the 
asymptotic behavior of problems in which the entries of at least one of the two matrices that specify 
the problem are chosen from a random distribution P. Surprisingly, this case has not been studied 
before using statistical methods despite the fact that the QAP was first proposed over 50 years ago 
01 . We find simple forms for Cmin and Cmax, the costs of the minimal and maximum solutions 
respectively. Notable features of our results are the symmetry of the results for Cmin and Cmax and 
the dependence on P only through its mean and standard deviation, independent of the details of 
P. After the asymptotic cost is determined for a given QAP problem, one can straightforwardly 
calculate the asymptotic cost of a QAP problem specified with a different random distribution P. 
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The quadratic assignment problem (QAP) is a combi- 
natorial optimization problem first introduced by Koop- 
mans and Beckmann 6] . It is NP-hard and is considered 
to be one of the most difficult problems to be be solved 
optimally. The problem was defined in the following con- 
text: A set of N facilities are to be located at N locations. 
The quantity of materials which flow between facilities i 
and j is Aij and the distance between locations i and j 
is Bij . The problem is to assign to each location a single 
facility so as to minimize (or maximize) the cost 
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where p{i) represents the location to which i is assigned. 

In addition to being important in its own right, the 
QAP includes such other combinatorial optimization 
problems as the traveling salesman problem and graph 
partitioning as special cases. There is an extensive lit- 
erature which addresses the QAP and is reviewed in 0- 
[isj . With the exception of specially constructed cases, 
optimal algorithms have solved only relatively small in- 
stances with < 36. Various heuristic approaches have 
been developed and applied to problems typically of size 
N ss 100 or less. 

Most work on the QAP has focused on solution tech- 
niques, bounds on optimal solutions, heuristics, and 
properties of problems with specially structured matri- 
ces. Previous work on asymptotic properties of ran- 
dom QAP instances has been limited to the case in 
which the elements of both matrices are drawn from ran- 
dom distributions. In this case it was shown by rigor- 
ous arguments [l9l - [23j that almost surely as — > oo, 
[Cmax — Cmin) I Cmax 0; thc minimum and maximum 
solutions approach thc solution obtained by a random 
permutation of p. 



Here we consider the properties of solutions to the 
QAP under the requirement that the elements of only 
one of the matrices need be drawn from random distribu- 
tion P. Our approach makes use of the replica approach 
of statistical mechanics. 

Without loss of generality, we will choose A as a ma- 
trix the elements of which are chosen from the random 
distribution P(Aij); the elements of B are arbitrary. We 
find that in the asymptotic limit in which the size of the 
problem N oo 

Cmin = iiaIJLbN^ -^Af{B)N^'^ (2) 
Cmax = iiaIIbN^ +^Af(B)N^'^. (3) 

Here Cmin and Cmax are the costs of the minimum and 
maximum solutions, respectively, and ^a and a a are the 
mean and standard deviations of the distribution P{A)] 
fiB the mean of the entries of B and 6 is a function of 
B and N. Our goal is to argue for the form of equations 
and ([3]). We do not attempt to determine the value 
of the functions f{B). 

It is useful to first consider the solution for which p=p* 
is a random permutation. Because the elements of A are 
assigned randomly and since p* (i) and p* (j) are random, 
each Bij in the sum is multiplied by a random value of 
Aij the average of which is ^a- Hence, the cost of a 
random permutation is 
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Bp*(i)p*{j) — f^Af^sN'^ 
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We now use the replica method of statistical mechan- 
ics to derive the form for Cmin and then derive the rela- 
tionship of Cmax to Cmin- Employing a Hamiltonian, H, 
defined as the QAP cost function our goal is to compute 
the partition function 
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and the free energy 
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where k and T are the Boltzmann constant and temper- 
ature respectively. Then, 



F{T = 0). 



Since the Hamiltonian includes a random matrix, A, we 
want to calculate the value of the free energy F averaged 
over the disorder specified by the probability distribution 
P{A). However, averaging the log of the partition func- 
tion is difficult. The replica method of statistical mechan- 
ics [13] was introduced to make calculation of this average 
possible. The replica method has been used not just on 
models of physical systems (such as spin glasses ;24-26|) 
but also on such combinatorial optimization problems 
as graph partitioning and the traveling salesman prob- 
lem [26 28]. The calculation of the average of the parti- 
tion function is simplified using a mathematical identity 
known as the replica trick, \n{x) = lim„_j.o(a;" — '^)/ri. 
Then equation (|6]) becomes 
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where 
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and = J P(A)Z'^dA denotes averaged over the 
disorder. Here each Hamiltonian represents a replica of 
the original system and the sum over {p"'} now denotes 
the sum over all permutations in all replicas. 

In order to achieve physically sensible results with / = 
F/N intensive, we require the following dependence on 
the mean and standard deviation of P(A) to scale as (see 

[2ii3): 
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with liA and a a independent of N . In Appendix A we 
then find that 
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We can make the following observations based on equa- 
tions (dll) and dll): 

• Consistent with equation the dependence of F 
on A is only through [ia and a a- 



• Vi OA = and/or T — )• oo, substituting equa- 
tion ([TTt in equation ([5]) yields F = ^AtJ-sN'^ 
which is the cost of the random solution, equation 
Q. This is reasonable because (i) physically for 
high temperatures we expect randomness and (ii) 
ii ga = 0, all entries in A are identical and all per- 
mutations yield the same costs. 



(7) We infer the form of F{T — 0) as follows: 



• In eauation (|lip ga appears in the combination 
ua/T. Thus, from equation ([8]) we see that in the 
T — >■ limit, only a term linear in a a can survive 
in F. 

This linear dependence on a a as well as on ^.a is 
consistent with the simple case in which all of the 
elements in A are scaled by a constant, z, in which 
case a A — >■ zga and ^a — > z^a- Clearly the opti- 
mal permutation is unchanged but the cost is also 
scaled by z. Thus, in this simple case, for any per- 
mutation (including the optimal one) the linear de- 
pendence on pLA and a a must hold. 

• Given that we obtained equation pT|) by expanding 
in 1/y/N, we expect the second term in the expres- 
sions for Cmin to be proportional N^/^ since the 
leading term is proportional to iV^. 

Given the considerations, above the only possible expres- 
sion for the second term in F is aAf{B)N^^^ where / is 
a function of B only. 

The form of Cmax follows directly as follows. Let 
Cmin(^, B) and Cma,x{A, B) denote the optimal minimum 
and maximum costs respectively of the QAP problem 
with matrices A and B and let — A denote a matrix with 
elements —Aij. Since Cinax(^, -B) = — Cinin(— ^, -B) and 
since ^-a = -^^^A and a -a = (^A, the form for Cmax in 
equation jS]) follows directly from the form for Cmin- 

If the entries of B are also drawn from a random dis- 
tribution, it is straightforward to show that 



Cmax = /iyiMsAf^ + CaAdBN 
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where c is a constant independent of A and B and (t_b is 
the standard deviation of the entries in B. 

Zdeborova et al. [2^ have conjectured that, in the 
large N limit, the minimal and maximal costs of parti- 
tioning random regular graphs into two equal sized sub- 
graphs are related by Cmax — = — Cmin where 
\E\ is the total number of edges in the random regular 
graph. Given that the graph partitioning problem can be 
represented as a QAP (see Appendix B) this relationship 
follows directly from equations ^ and ([3]) and supports 
this conjecture for the partitioning of random graphs into 
two subgraphs of any size. The more general relationship 

C^ax - ^^^^ = -H^Jif^ - C„„„ can be interpreted 
as an extension to weighted as well as unweighted graphs. 

We are not aware of a method to proceed further with 
the replica calculation for the case in which the B matrix 
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is not further specified, but if the matrix elements Bij 
can be represented as bibj (where bi, 621 ■• ■ can take 
on arbitrary values) the problem is tractable and can be 
taken further. This calculation is carried out in Appendix 
C and yields a result in the form of equation ([2]). 

To confirm our findings, we use the tabu search (TS) 
[soj heuristic to obtain approximate numerical solutions 
for a number of QAP instances. We employ matrices of 
the types described in detail in Appendix E. We use the 
notation "A matrix type"-"B matrix type" to specify a 
QAP instance. 

In Fig. [U we plot Cmin and Cmax versus aA for an 
instance of type Gaussian-Grid. As expected, the plots 
are linear in a a and the absolute values of Cmin and Cmax 
are equal for a given cr. This is consistent with equations 
© and ®. 

A stronger test is achieved by studying instances spec- 
ified by a matrix that represents a random graph of av- 
erage degree k. In this case, 
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{k) = y^k{N-l - k) 
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which represents a circle with origin at {{N — l)/2,0). 
In Fig. m^a) we plot Cmin, Crand, and Cmax versus k, 
< k < N — 1, for an instance of type Random-Grid. 



In order to illustrate the behavior of Cmin 
more detail, in Fig. [Hb), we plot 



and Cn 



rand 



^Cmin/max — Cmin/max C, 

The solid line is an ellipse of the form 



^theory ^ = ±(TA(fc)/(B)iV'/' 



(14) 



(15) 



where /(-B) is chosen to best fit of equation ([T5|) to the 
data. The fit is consistent with the theory, exhibiting 
both the expected linear dependence of the optimal costs 



on UA (k) and the symmetry represented by equations ([2|) 
and ([3]). In Fig. 21 we show similar plots for other 
varied QAP instances. Of particular interest are the plots 
for instances representing graph partitioning; the plots 
illustrate the confirmation of the conjecture of Ref. [2^ 
for both random and random regular graphs and also the 
validity for partitioning of graphs into unequal sized sets 
of vertices. 

We now study the dependence of AC on N. We treat 
instances in which the A matrix is random or random 
regular and consider different types of B matrix. To com- 
pare results for instances of different sizes, we define the 
normalized quantities ACn 
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AC 
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With this normalization we expect 



ACn 



CTA 
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In Fig. [2Ia) we plot ACnorm for various values of N 
for the Random-Grid instance. We confirm the 7V~^/^ 
dependence by plotting 



AC 



collapsed 



= ACn 



(18) 



in Fig. ISlJb). The collapse is consistent with equation 
(IT71) . Additional plots for other instance types are shown 
in Fig. [5] 

In summary, using the replica method of statistical 
analysis, we have found simple forms for the minimum 
and maximum costs of QAP problems in which at least 
one matrix is determined by a random distribution. 

We thank S. V. Buldyrcv for helpful discussions and 
the Defense Threat Reduction Agency (DTRA) for sup- 
port. 
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FIG. 1: For an TV = 100 QAP instance consisting of an A matrix with elements from a Gaussian distribution and a B matrix 
representing a two-dimensional grid, (from top to bottom), Cmax, Crand, and Cmin versus standard deviation a a. For a given 
Cmax and Cmin valucs are equidistant from Crand value. 
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FIG. 2: For an A'^ = 100 QAP instance consisting of an A matrix representing a random graph and B matrix representing a 
two-dimensional grid, (from top to bottom), (a) Cmax, Crand, and Cmin versus average degree k and (b) ACmax and ACminversus 
k. In this and all following figures, the upper and lower semi-circles are the ACmax and ACmin plots, respectively. The solid 
circular line represents the theoretical prediction. 




FIG. 3: For a QAP instance of type Random-Grid, (a) normalized AC versus normalized k for instance sizes A'^^100 (light 
gray); N= 225 (medium gray), and A'^=400 (black); (b) corresponding collapsed plots (see equation (|18p . 
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Appendix A: Integration over disorder 

In the following we retain only terms which do not van- 
ish in the N ^ oo limit. This is equivalent to retaining 
terms only to second order in Aij. Because we want to 
maintain the exponential form, we write 



where y.,j = J P(Aij) exp[^ J2l A„B. 



dAii - 1 



Expanding yij to second order in Aij we have: 



n N 



/I 
P{A) cxp[— ^ ^«j5p°(.)p°0)]rf^ 

= ^exp[ln / P(A)e-ETS'.j^°^'^'^>'°«p°«'dy4] 
{p°} 

= ^exp[ln]J / P(Ay)eT^^°^'J-^p°(')''°0)dAjj] 

{p°} ^J 



-1 B, 



{p°} 

^exp[^ln(l + t/,,)] (Al) 



r 1 " 

J P(A,,)[1 + —YAvB, 



2{kTy 
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where /^a2 is the second moment(around zero) of P. Expanding ln(l + yij) to the second order in yij and substituting 
equation (jA2p . we have 



ln(l + yy) - ^E^P"«P 



E/ -^p°(»)p° 



(i) 



2(fcT): 



where we have only retained terms to 0{Aij). Finally we have 



{p"} 



U) + 2(fcr)2(E^P°Wp°0'))^] 
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2{kT)'- 



■E(E^p= 



(«)p°0') 
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(A3) 



(A4) 



where we use the fact that j Spa (j)pa is independent Appendix B: Relationship to Graph Partitioning 
of permutation. 

The problem of partitioning a graph into two sub- 
graphs of size rN and {\—r)N with the minimum number 
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of edges between the two subgraphs can be represented where fj,2b is the second moment of the elements of the 
as a QAP as follows: One matrix, A, is the adjacency vector b, equation (IC1|) becomes 
matrix of the graph to be partitioned. The other matrix, 
B, the graph partitioning matrix, is the adjacency matrix 
for a bipartite graph in which edges are present between 
two sets of vertices; one set contains rN vertices and the 
second set contains (1 — r)N vertices. The QAP cost 
function is the cost of partitioning the graph represented 
by A 



Appendix C: Specific B Matrix 
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n N 



2{kT) 



N 



Here we treat the case in which the matrix elements 
Bij can be represented as hibj. We follow the spin-glass 
calculation of Ref . [1^ . 

Let 



N 



^'n ^ E ^-P[^(#)I E (E B,.,^,.,,)/Nf] (CI) 



{p°} 



a—1 



We can now use the Gaussian integral identity 



so 



(C2) 



Using 



N 



1 



'2ti 



E (E VwV(j))^ 

i.j — l a—1 

n N 

E (E^p°wVw)^ 

n N ^ 

nN'^fil^+ ^ 6pc(j)6p^(j))2, (C3) with dx { 2{kT)'i Y^'^dQ^p and find 

I 



ZL 



exp 



2b) 



al3 



-27r(fcr)2 



dQap exp [-ND[Qat3]] 



(C4) 



where 



D[Qo.0] 



Q^.-ln^expl(§r' 



N 



(a#;9) = l 1=1 



(C5) 



The form of equation ()C4p suggests that the integrals be evaluated with the method of steepest descent with the 
value of the integral determined by the maxim value of D. Assuming no replica symmetry breaking, at the maximum 
all values of Qap are equal [26j ; we denote this maximum value as Q and 



n N 



{P"} 



(C6) 
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Using 



n N 

(q//3) = 1 1=1 

N n n 

i—l Q a. 
N n 



EE(Vw)'-^^M26 (C7) 



to uncouple the replicas, we have 



N n 



D[Q] = -(^)'("(^ - + 2^A^2.Q) - In E exp(g)'^ E(E (C8) 



i— 1 a 



In the spin glass calculation the sum over i can be carried out here and the factor of drops out. Since we cannot 
perform the sum over i here, the 1/A'' factor is carried through to the end of the calculation as a normalization factor. 
Using the Gaussian integral identity again we can write 

^ N n 

E[Q] ^ ln^exp[(^)2-^(^6p=(,)fQ/7V] 

= 1- E n / --p[-| + §Vq7a^E v«-.]4= (C9) 



Now considering just the terms in the exponent dependent on fepQ^j), we find 

~ n ~ ~ n 

E exp[^VQ7^E^f°w^^] = n(Ee^p[^^^Vw^»]) = (Eexp[gyo7A^E^pw^'])" (CIO) 

This is the key step which makes the n dependence explicit and allows non-ambiguous analytic continuation of n — 0. 
Keeping only terms to first order in n, and using (l/27r) e^^ /'^dz = 1 we can then write 

E[Q] = InH / 6-^^^/2(5] exp[^yQ7iVV0^»])"7=) 

i P 

= InH / e-?/2(i + ^in^exp[^V^&pW^.])7=) 
= ln(l + nn / e-?/2ln5]exp[^v^6pW^.]^) 

= nf] / e---/2in^exp[^v^6.W^^]7= (Cll) 
For small T, as shown in Appendix [Dl 



m^n^VO/N^ib), (C12) 



where depends only on the values of bi normalized by y/N. We can then write 



D[Q] = -(^)'(«(" - + 2n/.2fcg) - (C13) 
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To find the Q which maximizes Z' we differentiate equation (jC13|) with respect to Q and solve for Q. In the limit 



n 



we find 



kT 

Q = m-^^- (C14) 
Substituting for D\Q] and evaluating equation (jC4p using the method of steepest descent we find, to first order in 

= eM-N{-^f{-^4bn + n{n-l)Q^ + 2n^i2bQ)-n^^■^{b)] 

= eMnN{{^f{{^l^, - ~ ^^^{h))] (C15) 



and thus from equation ()C2p 

F{T^O) = - kT \im -(Z^ ~1) 

n-i-O n 

= -kT[^iA^iBN^ + (^)'^((m26 - Qf - ^VQ^m (Ci6) 

Using equation (|C14p . in the limit T — > 0, and defining f{b) = ^fi2b'^{b) 

F(T = 0) = fiA^isN^-aANfib) 

= ^lA^iBN^-aAN^/^f{b) (C17) 



Appendix D: E[Q] calculation 

Here we discuss the evaluation of We want to evaluate multiple integrals of the form 

E[Q]=ll / e--/2ln^cxp[^VQ7^6pW-.]^- (Dl) 
To see how the evaluation would proceed, consider the case of N=2. Then 

where a ^ 1/T . For small T, for different regions of integration, one exponential in the sum of exponentials dominates. 
Specifically, assuming without loss of generality that bi > 62, 

l^^^a{b^.^+b2Z2 ^ ^a{b2z,+b^z2^ _ abiZi + \n[e''^^'^] (Zi > Z2) 

- 062^1 +ln[e''''i^=] {zi<Z2). (D3) 

Then 

E[Q] ^ J e-^'/^abiZi + \n[e''''''']dzl + J'' 062^1 + In^''^^^]]-^)-^. (D4) 

These integrals can all be solved exactly, with intermediate results in terms of the error function, erf{x) = 
e~^^dz, and we find that for this case of iV = 2,E[Q] = a{bi — b2)/\/n. This approach can be extended 
to any N and the result is a times a linear combination of the constants bi . We can then write 

E[q]^n^^^ib), (D5) 



where ^(6) is the linear combination of the bi divided by which normalizes the expression. 

Appendix E: Matrix Types • Uniform - the matrix elements are chosen from a 



We employ matrices of the following types: 
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uniform distribution on the interval [0, 100]. 

• Gaussian - matrix elements are chosen from a Gaus- 
sian distribution with zero mean and standard de- 
viation cr. 

• Half-Gaussian - matrix elements are chosen from a 
Gaussian distribution as above but only elements 
with value greater or equal to zero are used. 

• Random (graph) - the matrix is the adjacency ma- 
trix of a random graph with edges present with 
probability p. The average degree of the graph is 
k = pN. 

• Random Regular (graph) - the matrix is the adja- 



cency matrix of a random regular graph for which 
all vertices are degree k. 

• Grid - the matrix elements are the Euclidean dis- 
tances between points in a two-dimensional square 
grid. The distances between adjacent points along 
the X and y axes are 100. 

• Graph Partitioning - the matrix is the graph parti- 
tioning matrix described in Appendix [Bl 

All matrices are symmetrical with zero diagonal. For the 
Random and Random Regular matrices that represent 
graphs, we study cases of the graph degree k ranging 
from to — 1. 




FIG. 4: For various A*' — 100 QAP instance, ACmax and ACminversus k. The solid circular line represents the theoretical 
prediction. 
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FIG. 5: (a),(c),(e) Normalized AC versus normalized k for instance sizes Ar=100 (light gray); N= 225 for (a )and 200 for (c) 
and (e) (medium gray); and A'^=400 (black). The right hand column contains the corresponding collapsed plots. Panels (a) 
and (c) also appear in the main paper. 



